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Abstract Past and modern large-scale ice sheet loss results in geographically variable sea level 
changes. At present, in Hudson Bay, Canada, sea level is decreasing due to glacial isostatic adjustment, 
which represents a departure from the globally averaged sea level rise. However, there are large 
uncertainties in future sea level trends with further polar ice sheet retreat in the coming centuries. Sea 
level changes affect ocean tides considerably because tides are highly sensitive to changes in bathymetry. 
Here, we present multi-century sea level projections associated with a suite of past and future ice loss 
scenarios and consider the impact of these changes on ocean tides using an established tidal model. 
Modern tides in Hudson Bay are poorly resolved due to large uncertainties in bathymetry. To establish an 
initial condition for our simulations, we constrain bathymetry in the bay using tide observations. Due to 
gravitational, Earth rotational and deformational effects, Greenland ice loss will produce a small sea level 
fall in the bay, while Antarctic ice loss will produce a larger than average sea level rise. Our results show 
that the response of the Antarctic ice sheet to climate change strongly impacts the magnitude and sign of 
future sea level and tidal amplitude changes in the region, with the largest changes predicted in Hudson 
Strait and Foxe Basin. We emphasize that further constraints on bathymetry and accurate projections 

of sea level and tides in Hudson Bay are imperative for assessing the associated impacts on coastal 
communities and ecosystems. 


Plain Language Summary Hudson Bay is a shallow bay in northern Canada surrounded 

by coastal communities and ecosystems that are vulnerable to future sea level change. The bay was ice 
covered 21,000 years ago, and sea level is currently falling there due to ongoing land uplift since the 

ice retreated. It is unclear if this trend will continue as the Greenland and Antarctic ice sheets melt, 
contributing to spatially variable sea level changes. Sea level changes also impact ocean tides due to their 
sensitivity to water depth. We model future sea level and tide changes in the Hudson Bay region associated 
with land uplift and projections of Greenland and Antarctic ice loss over the next 500 years. With less 
Antarctic ice loss, sea level continues to fall and tidal amplitudes decrease. Conversely, high-end projected 
Antarctic ice loss could increase water depths in Hudson Bay by up to 14.5 m by 2500 and tidal amplitude 
changes could exceed 1 m, increasing in much of Hudson Bay while decreasing in Hudson Strait. A better 
understanding of the response of Antarctica to climate change will improve projections of sea level and 
tide changes in the Arctic and the associated societal and environmental impacts. 


1. Introduction 


Global mean sea level rise has accelerated in recent decades (Chen et al., 2017; Hay et al., 2015; Nerem 
et al., 2018), and this acceleration is expected to continue as global temperatures rise, resulting in thermal 
expansion of the oceans and melting of mountain glaciers and the polar ice sheets (Church et al., 2013; 
Clark et al., 2015; Slangen et al., 2016). Of these contributors, the Greenland and Antarctic ice sheets hold 
the largest potential to increase global mean sea level (WRCP Global Sea Level Budget Group, 2018), and are 
expected to be the dominant contributors to global mean sea level rise on multi-century timescales (Church 
et al., 2013; Golledge et al., 2015). Regional sea level changes associated with ice loss can differ substan- 
tially from the global mean, due to gravitational, Earth deformational, and rotational effects (e.g., Clark & 
Lingle, 1977; Gomez et al., 2010; Mitrovica et al., 2001). Furthermore, past (ice age) ice cover changes lead 
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Figure 1. Model domain used in this study, showing (a) Bathymetry and topography of the Hudson Bay Complex (data from the composite dataset described 
in Section 3) defined here as Hudson Bay, James Bay, Foxe Basin, and Hudson Strait. (b) Present day semi-diurnal (M2) tidal amplitudes (color) and phases 
(contoured in white lines at 1/8 M, period), taken from the TPXO8 tide database (http://volkov.oce.orst.edu/tides/tpxo8_atlas.html). 


to ongoing glacial isostatic adjustment (GIA) that contributes to regional variability in sea level change. 
GIA effects can be considerable in regions such as the Hudson Bay Complex (HBC, defined here as Hudson 
Bay, James Bay, Foxe Basin, and Hudson Strait, Figure 1a) in northern Canada that were once ice covered 
during the Last Glacial Maximum. In the HBC, sea level is presently falling due to GIA, but it is unclear 
if this trend will continue in the coming centuries due to the large uncertainties in the future evolution of 
the polar ice sheets, in particular the Antarctic ice sheet (e.g., Church et al., 2013; DeConto & Pollard, 2016; 
Golledge et al., 2015). 


Here, we define sea level as the height of the sea surface relative to the solid Earth surface, or equivalently, 
ocean bathymetry. Global and regional patterns of sea level change have been shown to have a strong in- 
fluence on ocean tides (e.g., Padman et al. 2018; Wilmes et al., 2017), because tides travel as shallow water 
waves, which are sensitive to ocean bathymetry. Changes in ocean bathymetry (caused by changes in sea 
surface height and/or the shape of the sea floor) can lead to changes in tidal wave propagation speed, en- 
ergy dissipation, and amplitudes (Green, 2010). Variations in sea level and the associated changes in ocean 
tides can, in turn, affect other coastal processes and ecosystems such as the development of salt marshes, 
coastal erosion, frequency and magnitude of flood events, and saltwater intrusion of surface water (Craft 
et al., 2009; Kirwan et al., 2010; Kirwan & Guntenspergen, 2010; Kirwan & Megonigal, 2013; Nicholls & 
Cazenave, 2010; Ross et al., 2017). Finally, Arctic and sub-Arctic indigenous communities surrounding the 
HBC are among the most vulnerable to climate change (IPCC, 2014), with their livelihoods and food sources 
critically impacted by coastal changes (e.g., Lemmen et al., 2016; Tsuji et al., 2009, 2016). 


Projections of future sea level change in the HBC due to climate change have been relatively limited. 
Gough (1998) suggested that sea level in Hudson Bay could begin to increase in the current century as the 
global mean sea level rise associated with climate change outpaces the sea level fall due to glacial isostatic 
uplift. However, Gough's simple model did not capture the spatial variability in sea level change that arises 
from melting of current ice sheets. More recently, Tsuji et al. (2009, 2016) and Lemmen et al. (2016) calcu- 
lated the spatially variable sea level change and suggested that glacial isostatic uplift would persist in parts 
of the HBC in the coming centuries but did not consider a wide range of possible climate-driven sea level 
changes. Recent work has shown that even under the same climate forcing scenario, Representative Con- 
centration Pathway (RCP) 8.5, projections of Antarctica's contribution to global mean sea level rise range 
from 2.25 to 15.65 m by 2500 (DeConto & Pollard, 2016; Golledge et al., 2015, 2019). Given this large uncer- 
tainty, it is unclear what the dominant contributor to sea level change in the HBC will be, or even what the 
sign of sea level change will be. Furthermore, due to the relatively shallow bathymetry of the HBC, and the 
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low-lying topography that surrounds the bay (Figure 1a), a small rise or fall in sea level will result in a large 
shoreline migration, amplifying the HBC's sensitivity to uncertainty in future sea level change. 


Previous regional studies have investigated future sea level and tide changes on the European Shelf (e.g., 
Idier et al., 2017; Pelling & Green, 2014; Pickering et al., 2012), the Patagonian Shelf (Carless et al., 2016), 
the Gulf of Maine (Pelling & Green, 2013), the Gulf of Mexico (Passeri et al., 2016), and the Bohai Sea in 
China (Pelling et al., 2013), but studies of the HBC are lacking. The HBC is a particularly interesting region 
with regards to ocean tides. The tidal range of Ungava Bay in Hudson Strait (16.8 + 0.2 m, Drinkwater, 1986) 
rivals that of the Bay of Fundy (17.0 + 0.2 m, O'Reilly et al., 2005) as the region with the largest tidal range 
in the world. Furthermore, Egbert and Ray (2001) found that more tidal energy is dissipated in the HBC 
alone than in any other region in the world. This is because the HBC is near resonant at the semi-diurnal 
tidal forcing frequency (Arbic et al., 2007; Webb, 2014), producing large semi-diurnal tides (exceeding 4 m 
amplitudes in some areas, Figure 1b) that are further amplified due to the HBC's location on a continental 
shelf (Clarke & Battisti, 1981). The resonant properties of the HBC and its shallow continental shelf sea 
character can strongly influence deep ocean tides in the connecting Atlantic Ocean (Arbic & Garrett, 2010; 
Arbic et al., 2009), so any large-scale change in the HBC can have consequences far away from the region. 
Previous work on paleotides in the region indicated a strong sensitivity of ocean tides to changes in sea level 
(Arbic et al., 2004; Egbert et al., 2004; Uehara et al., 2006). However, modern tides in the HBC are not well 
resolved due to limited bathymetry data that can be used to constrain models (Egbert & Ray, 2001). Project- 
ing future tides therefore requires, as a first step, an appraisal of the performance of available bathymetry 
models in reproducing modern tides in the HBC. 


Wilmes et al. (2017) showed that a future complete collapse of the Greenland and West Antarctic ice sheets 
would strongly impact tides globally, and that simulations that take into account spatial variability in sea 
level change differ substantially from those that assume the global average equivalent rise. In particular, in 
Hudson Bay, the spatially variable sea level change associated with Greenland ice loss is opposite in sign to 
the uniform global average sea level rise scenario, which, in turn, results in large differences in predicted 
tidal changes. Furthermore, they show that spatially variable changes in sea level, and in turn tidal ampli- 
tudes in much of the HBC are positive under West Antarctic ice sheet collapse and negative under Green- 
land ice sheet collapse. 


As a complete collapse of the Greenland and West Antarctic ice sheets represents a long term, upper bound 
on their contribution to sea level change, a question arises: How will future sea level and tides in the HBC 
evolve under combined, transient melting of the ice sheets in the coming centuries? Greenland's contri- 
bution to sea level change, which is not expected to exceed 20 cm of global mean sea level equivalent by 
2100 (Church et al., 2013; IPCC, 2019), is predominantly driven by surface mass balance, followed by ice 
discharge via icebergs (Church et al., 2013; IPCC, 2019). The HBC is relatively insensitive to Greenland ice 
loss due to its close proximity to the net-zero sea level change line associated with gravitational effects that 
arise from Greenland melting. Conversely, the dynamic response is expected to dominate sea level changes 
due to Antarctic ice loss, but the future contribution from Antarctica remains highly uncertain due to the 
lack of observations to constrain models and an incomplete description of ice sheet-ocean-atmosphere in- 
teractions in current models (IPCC, 2019; Kopp et al., 2017). Projections range from tens of centimeters to 
more than 10 m global mean sea level equivalent over multi-century timescales (e.g., Golledge et al., 2019; 
Pollard et al., 2017). In contrast to Greenland melting, gravitational effects due to Antarctic ice loss will lead 
to a larger-than-global-average sea level rise in the HBC. Here, we investigate to the opposing influence 
of gravitational effects and the large uncertainty in Antarctic projections on sea level and tide changes in 
the HBC. We consider a single contribution from the Greenland ice sheet and low- and high-end possible 
contributions to sea level change from the Antarctic ice sheet under RCP8.5 climate warming. On decadal 
timescales, sea level changes will likely be dominated by GIA, but on centennial timescales, the contribu- 
tion from ice sheets could exceed GIA, in particular under high-end warming scenarios. This transition 
from a regime of sea level fall to one of rise, and the associated changes in tides, is important to consider in 
assessing coastal climate change impacts and develop adaption strategies. 


The study is structured as follows. In Section 2, we describe the ice loss and sea level projections and the 
tide model setup. Due to the sensitivity of tides to bathymetry, we discuss constraints on modern bathym- 
etry in the HBC in Section 3 and present a bathymetry model that improves modern nonassimilative tide 
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simulations. In Section 4, we present predictions of future sea level changes associated with future retreat 
of the polar ice sheets and ongoing GIA and the associated tidal changes in the HBC. We conclude with a 
discussion in Section 5. 


2. Methods 
2.1. Sea Level and Ice Sheet Modeling 


We adopt the sea level theory from Kendall et al. (2005) and Gomez et al. (2010) to predict the spatially 
variable sea level changes resulting from past and future ice cover variations using a gravitationally self-con- 
sistent sea level model that incorporates viscoelastic deformation of the solid earth due to surface ice and 
ocean mass loading, migrating shorelines, and Earth rotational effects. Elastic and density structure of the 
solid Earth in all simulations is provided by the Preliminary Reference Earth Model (Dziewonski & An- 
derson, 1981) and viscosity structure is given by the VM2 Earth model (Peltier, 2004) unless otherwise 
specified. The modern topography that serves as input to the sea level model is taken from GEBCO 2014 
(Weatherall et al., 2015) globally, with the composite regional bathymetry described below over the HBC. 
All calculations are performed up to spherical harmonic degree 512. 


We simulate the future evolution of the polar ice sheets in Greenland and Antarctica with dynamic ice sheet 
models. We focus on projections adopting RCP 8.5 (Collins et al., 2013), representing the business-as-usual 
emissions scenario as it serves as the upper bound of expected global temperature increase, and hence 
contribution of the ice sheets to sea level change in the HBC. Furthermore, by considering a single RCP 
scenario, we highlight the variability that emerges among ice loss and sea level projections that are guided 
by the same emissions framework. 


Greenland and Antarctic ice thickness projections are generated using the Parallel Ice Sheet Model version 
0.7.3, a hybrid ice model with shallow shelf-shallow ice approximations for floating and grounded ice (Win- 
kelmann et al., 2011). Following the procedure outlined in Golledge et al. (2019), future climate forcing 
on the ice sheets is driven by the Coupled Model Inter-comparison Project Phase 5 multi-model ensemble 
mean outputs. The ice sheet model simulations are performed on a high-resolution (2.5 km for Greenland 
and 5 km for Antarctica) polar stereographic grid. We apply a Gaussian smoothing filter on the changes in 
total (floating and grounded) ice thicknesses, then interpolate ice thickness predictions onto a lower resolu- 
tion Gauss-Legendre 512 x 1024 global grid, to serve as input to the sea level model. We consider the Antarc- 
tic projection from Golledge et al. (2019) as the low-end projection of future Antarctic ice loss under RCP8.5 
(henceforth the “low-end” scenario). We also consider an alternative simulation for the Antarctic ice sheet 
that predicts much greater ice loss under RCP8.5. This high-end projection of Antarctic ice thickness chang- 
es (hereafter referred to as the “high-end” scenario) is provided by a simulation in Pollard et al. (2017) using 
a coupled ice sheet—sea level model that includes ice cliff and hydrofracturing physics and a hybrid combi- 
nation of shallow ice-shallow shelf approximations for ice dynamics. Ice thickness variations are provided 
by the ice model at 20 km resolution and then linearly interpolated onto the global 512 x 1024 grid. 


To fully quantify the sea level changes in Hudson Bay, we also incorporate the contribution of ongoing GIA 
associated with past ice cover changes over the last deglaciation, which, in our simulations, spans the time 
period of 21,000 years ago to the year 2000 (henceforth the modern). We model ice loading changes over 
the last deglaciation by adopting the ICE-5G ice history with the VM2 Earth model (Peltier, 2004) and the 
ICE-6G ice history with the VM5a Earth model (Argus et al., 2014; Peltier et al., 2015). Each simulation was 
run from 21,000 years ago to 500 years in the future, and the predictions of change from the modern to the 
year 2500 were added to the sea level changes associated with future ice sheet retreat described above (see 
Figure S1 for a comparison of ICE-5G and ICE-6G). 


The total sea level change is calculated as the sum of the contributions from Greenland, Antarctica, and 
GIA. Note that our projections do not include oceanographic effects. In the current century, these are ex- 
pected to contribute a relatively uniform sea level rise in the bay under future climate warming (Jevrejeva 
et al., 2016), but the magnitude, timing and spatial variability of these changes is challenging to simulate in 
the HBC on multi-century timescales. While these effects would alter the timing of the change from a sea 
level fall to a sea level rise, they would likely not offset the overall effect of transitioning from a sea level fall 
to a sea level rise. In addition, the size of the contribution is well within the range of Antarctic contributions 
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considered here. Shorter timescale simulations should include the influence of ocean dynamics on sea level 
and tides. 


2.2. Tide Modeling 


We use the Oregon State University Tidal Inversion Software (OTIS; Egbert et al., 1994, 2004) to model 
present-day and future tides in Hudson Bay. OTIS solves the linearized shallow water equations forced by 
the tidal potential and with energy dissipation from bed friction and parameterized tidal conversion (e.g., 
Green & Nycander, 2013; see Wilmes et al., 2017 for details on the model setup). 


The HBC tide model was run at 1/30° x 1/30° horizontal resolution and simulated the tidal constituents 
M>, S2, and K;. The extent of the model domain is shown in Figure 1. Dissipation of tidal energy was cal- 
culated according to Egbert and Ray (2001; see also Wilmes et al., 2017 for a detailed description). For the 
present-day simulations, TPXO8 (Egbert & Erofeeva, 2002; see http://volkov.oce.orst.edu/tides/tpxo8_atlas. 
html for the latest version) elevations were used as forcing at the open boundaries. TPXO8 is an observation 
constrained tide model solution that assimilates tide gauge observations and satellite altimetry to produce 
a highly accurate estimate of the tides (~2 cm RMS error compared to unassimilated tide gauge data; see 
Egbert & Erofeeva, 2002, and Stammer et al., 2014, for details). 


For the future tide simulations, sea level change predictions described in Section 2.1 were added to the pres- 
ent-day bathymetries for both HBC and the global setup. Simulations were performed for 2100, 2300, and 
2500 under low- and high-end Antarctic contributions. Elevation changes in response to each of the future 
sea level change scenarios simulated by a near-global (86°S-89°N) tidal model (see next paragraph) were 
added to the boundary conditions used for the present-day simulations, thereby ensuring a high degree of 
accuracy for high-resolution HBC future simulations. We carried out two sets of runs for each time slice: 
one where flooding of low-lying land was permitted (denoted “flood”) and a second one where flooding was 
not permitted (“no-flood”). 


The near-global tide model simulations were also run with OTIS (see Wilmes et al., 2017 for details on lin- 
earized shallow wave equations, and internal drag and self-attraction and loading parameterizations). The 
underlying bathymetry was RTopo-2 (Schaeffer et al., 2016) and the model was run at 1/8° x 1/8° horizontal 
resolution for 21 days to simulate M,, S,, and K;. The M, amplitude errors when compared to TXPO8 are 
5.6 cm globally and 2.8 cm in the deep ocean (see Table S1). Integrated M, tidal dissipation for the global 
and deep ocean is within 6% and 7%, respectively, of values computed for TPXO8, despite the very different 
resolutions in the shelf seas (TXPO8 has a resolution of 1/30 degree in most shelf seas; see also Egbert & 
Ray, 2001). Simulations using the same sea level scenarios as for the regional model were carried out with 
the near-global model and used as boundary conditions for the regional model simulations. Land-ocean 
boundaries in this case were kept constant (i.e., land was not permitted to flood) due to the low resolution 
of the model. 


3. Constraining Modern Bathymetry 


Numerical prediction of tides in near-resonant basins, such as Hudson Bay, are strongly sensitive to the 
adopted ocean bathymetry (Egbert & Ray, 2001; Green, 2010; Padman et al., 2018). Therefore, in order 
to project future tide changes, we require an accurate estimate of modern bathymetry. Bathymetry in the 
HBC is poorly constrained due to the sparsity of available data, related to the HBC's shallow depths, spatial 
extent, and limited commercial shipping routes. Furthermore, there is substantial disagreement between 
global bathymetry datasets in the HBC (Figure S2); differences between datasets reach up to 200 m in some 
areas of the complex. As a result, tide models struggle to accurately capture the tides in the region. 


Given the lack of available data to distinguish between these bathymetry models, we seek to simultaneously 
find a bathymetry to use for (future) projections and to improve non-assimilative modeling of the modern 
tides in HBC. To do this, we simulate modern tides in the HBC using the configuration described in Section 
2.2 with a suite of different bathymetry models GEBCO 2008 (see The GEBCO_2008 Grid, version 20100927, 
www.gebco.net), GEBCO 2014 (Weatherall et al., 2015), ETOPO1 (Amante & Eakins, 2009), and SRTM30_ 
Plus (Becker et al., 2009) and evaluate their performance, focusing on Hudson Bay and Strait. 
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Figure 2. Amplitude differences between simulated M, amplitudes from several bathymetry datasets and M, amplitudes from the TPXO8 tidal solution. (a) 
ETOPO1 (Amante & Eakins, 2009), (b) GEBCO 2008 (http://www.gebco.net), (c) SRTM30_Plus (Becker et al., 2009), (d) GEBCO 2014 (http://www.gebco.net), 
(e) our new, composite bathymetry with SRTM30_Plus in Hudson Bay and GEBCO 2008 elsewhere. The colored circles show the amplitude differences for the 
tide-gauge locations (see also Supporting Information S1 and Table S4). 
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Figures 2a-2d show the difference between the simulated M, tidal amplitudes using the four bathymetry 
datasets considered (ETOPO1, GEBCO 2008; SRTM30_Plus, and GEBCO 2014) and the TPXO8 amplitudes. 
We calculate the root mean square complex elevation errors (RMSEs) against 14 tide gauges (see Supporting 
Information S1 for details and Figure 2 for their locations), TPXO8 and TPXO9 (see https://www.tpxo.net/ 
global/tpxo9-atlas for details and the latest version) to evaluate the performance of bathymetry models in 
reproducing modern tides over the global domain (see Table 1 for a quantification of the errors). Amplitude 
and phase RMSEs are shown in Tables S2 and S3 for all ocean depths (total RMSE) and deep oceans regions 
(water depths > 500 m, deep RMSE). 


The simulation using ETOPO1 produced the largest elevation and amplitude errors of all the simulations 
(Figures 2a and Table 1), with amplitudes differing from TPXO8 by more than 1 m in Hudson Strait and 
James Bay. Both GEBCO 2008 (Figure 2b) and GEBCO 2014 (Figure 2d) perform well in Hudson Strait, 
each with mean amplitude differences of 21 and 15 cm in Hudson Strait, respectively. In Hudson Bay, 
SRTM30_Plus (Figure 2c) matches more closely than the GEBCO simulations, differing from TPX08 with a 
basin amplitude RMSE of 7 cm. The largest amplitude errors in all simulations are in Hudson Strait, where 
there are also the largest differences in depths between datasets (Figure S2) and the largest tidal amplitudes 
(Figure 1b). 


To improve the accuracy of our nonassimilative tidal simulations, we develop a composite bathymetry (Fig- 
ure 2e) based on how well each existing bathymetry dataset performed regionally in the comparison to 
TPXO8 in the HBC. Our composite bathymetry merges, using linear weighting, SRTM30_Plus bathymetry 
in Hudson Bay and GEBCO 2008 bathymetry elsewhere in the HBC. The amplitudes of the tide simulation 
adopting the composite bathymetry differ from TPX08 in Hudson Bay on average by 6 cm. The largest am- 
plitude errors in James Bay and Foxe Basin are less than +25 cm, with mean errors of 5 and 13 cm, respec- 
tively. Our composite bathymetry improves the amplitude fit in Hudson Strait, with an average amplitude 
difference of 20 cm when compared to TPXO8. Furthermore, we find that our composite, regional bathym- 
etry substantially reduces elevation errors both in comparison to tide gauges, TXPO8 and TXPO9 (Table 1) 
and also both total and deep RMSE for amplitudes and phases (Table S3) for all simulated constituents. 
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Consequently, we adopt this composite, modern bathymetry in the future 
sea level and tide projections shown in the sections to follow. 


ETOPO1 
HB 
HS/FB/CS 
SRTM 
HB 
HS/FB/CS 
GEBCO 2008 
HB 
HB/FB/CS 
GEBCO 2014 
HB 
HB/FB/CS 
COMPOSITE 
HB 
HS/FB/CS 


70.7/36.7/38.3 
39.2 
127.8 
50.2/33.9/37.8 
56.6 
38.5 
54.3/24.2/25.3 
38.3 
83.3 
47.8/29.8/30.9 
25.8 
87.4 
19.3/14.4/15.9 
15.0 
27.0 


S2 Kı 
25/10.8/11.8 4.7/3.2/4.0 4. Results 
12:7 3.6 
4.1. Sea Level Change 

47.1 6.5 
22.2/10.3/11.6  3.9/2.8/3.3 Figure 3 shows the projected contributions to sea level change in the 
Dis AA HBC associated with GIA and melting from the Greenland and Antarctic 
i i ice sheet for the three time slices 2100, 2300, and 2500. The contribution 
Aad 28 from GIA (Figures 3a—3c) shows a sea level fall across the HBC. James 
17.4/8.4/9.2 3.0/3.3/4.1 Bay and Churchill, Manitoba regions experience the maximum predicted 
10.1 32 sea level fall (see Figure 1a), with a magnitude of 1.2 m in both regions 
30.4 25 by 2100, and 5.9 and 5.7 m respectively, by 2500. In Foxe Basin, sea lev- 
16.6/9.2/10.1 3.6/2.8/3.3 €! fall associated with GIA reaches 0.9 m by 2100, and 4.3 m by 2500. 
BG oy The sea level fall is driven by the larger thickness of the Laurentide Ice 
i ; Sheet during the Last Glacial Maximum in these regions and the later 
ey a2 termination of deglaciation. Hudson Strait sees a more moderate fall of 
11.8/4.7/5.7 2.9/2.6/3.6 0.3 m by 2100 and 1.5 m by 2500. In contrast, to the east of the HBC in 
6.6 3.1 Davis Strait region, sea level rises over the 500-year GIA simulation (red 
21.2 2.4 regions in Figures 3a—3c) due to the subsidence of peripheral bulges that 


formed around the ice sheets at the Last Glacial Maximum in our ICE-5G 


Note. Shown are complex elevation errors in cm between the tide model 


simulations using the respective bathymetry and tide gauges/TPXO8/ 
TPXO9 for M2, Sz, and K;. The best fit has been indicated as bold for 
each constituent and evaluation method. For each bathymetry we also 
show errors against the tide gauges located in Hudson Bay (HB) only and 
the tide gauges located in Hudson Strait (HS), Frobisher Bay (FB), and 


Cumberland Sound (CS). 
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simulation. 


Next, we consider the contribution from melting of the polar ice sheets 
(Figures 3d-31). There are a number of physical effects contributing to 
the differences in sea level change in the HBC associated with Greenland 
(Figures 3d-3f) and Antarctic (Figures 3g-31) ice loss. An ice sheet exerts 
a gravitational attraction on the surrounding water, and this attraction 
weakens as the ice sheet loses mass, resulting in a local sea level fall with- 
in ~2,000 km of the ice loss (Woodward, 1888). At greater distances from 
the source of ice loss, sea level rises more than the global average (Tamisiea & Mitrovica, 2011). Viscoelas- 
tic uplift of the solid Earth in and around the region of ice loss further accentuates the local sea level fall, 
and deformation of the Earth due to loading of the oceans with water and changes in the planet's rotation 
associated with surface mass redistribution add geographic variability to the sea level change farther away 
from the ice sheet. In our simulations, the Greenland ice sheet is projected to contribute 0.10 and 0.55 m 
of global averaged sea level rise by 2100 and 2500, respectively, but the HBC is within the zone of gravi- 
ty-driven sea level fall for the Greenland ice sheet. Therefore, Greenland melting is projected to contribute a 
low-amplitude drop in sea level across the HBC that decreases in amplitude with distance from Greenland. 
The smallest sea level decrease occurs in the south-eastern part of the HBC (Figures 3d-3f and S3a) and the 
largest sea level fall can be seen in northern Foxe Basin (up to 0.7 m by 2500) and the northern Labrador 
Sea. In Hudson Strait, the projected sea level decrease reaches 0.5 m by 2500. Central Hudson Bay would ex- 
perience up to a 0.2 m decrease in sea level by 2500. In James Bay, sea level remains unaltered by 2100, with 
a shift to a sea level rise by 2500 that does not exceed 0.10 m. The sea level fall due to Greenland melting is 
an order of magnitude smaller than that caused by GIA. Sea level in the HBC is hence relatively insensitive 
to Greenland ice loss due to the fingerprint pattern and the proximity of the southern and western coasts of 
the HBC to the region of net-zero sea level change. 


Melting from the Antarctic ice sheet, on the other hand, would lead to a sea level rise in the HBC that is, 
approximately 1.2 times greater than the global mean sea level equivalent rise associated with the melting, 
with a relatively uniform rise across the HBC (Figure S3b). Figures 3g—31 show projected sea level changes 
in response to the low-end (Figures 3g-3i) and high-end (Figures 3j-gl) responses of the Antarctic ice sheet 
to the RCP8.5 scenario. Under the low-end simulation, Antarctic melting would lead to a maximum sea 
level rise in the HBC of 0.1 m by 2100 and 2.1 m by 2500. The high-end scenario (j-1) predicts an order of 
magnitude larger sea level rise of up to 0.8 m by 2100 and 14.5 m by 2500 across the HBC. 
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Figure 3. Contributions to sea level change in the Hudson Bay Complex from glacial isostatic adjustment (GIA) and future melting of the polar ice sheets 
under RCP8.5 at 2100, 2300, and 2500, relative to 2000. Blue shading corresponds to a sea level fall and red shading corresponds to a sea level rise. (a-c) 
Contribution to sea level change from past ice-ocean loading changes (i.e., GIA) over the last deglaciation associated with the ICE-5G ice history (Peltier, 2004). 
(d-f) Contribution to sea level change from the Greenland ice sheet (Golledge et al., 2019). (g-1) Contribution to sea level change from the Antarctic ice sheet 
under from (g-i) low-end (Golledge et al., 2019) and (j-1) high-end (Pollard et al., 2017) projections. Note the different color scales, the color bar associated with 
each panel is on the right side of each plot. 
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In Figure 4, we consider the sum of all the contributions shown in Figure 3 (see also Figure S4 for the total 
sea level change adopting ICE-6G). Under the low-end Antarctic ice loss scenario (Figures 4a-4c), the sea 
level changes at 2100, 2300, and 2500 are dominated by GIA effects. The largest sea level decreases are pre- 
dicted in the regions where the signal from GIA is greatest (Figures 3a-3c). At 2100 (Figure 4a), the magni- 
tude of sea level fall ranges from 1.1 m in Churchill to 1.2 m in James Bay. In Hudson Strait, the amplitude 
of sea level fall is damped and on the order of 0.4 m. East of Hudson Strait in the Davis Strait, a slight rise 
is seen, where Antarctic ice loss and GIA are both contributing a sea level rise. By 2300 (Figure 4b), the sea 
level fall ranges from 1 m in central Hudson Bay up to 2.75 m in James Bay. Ungava Bay, at the entrance 
of Hudson Strait, is the only area in the HBC that experiences a sea level rise, reaching up to 0.4 m. Thus, 
although Greenland ice loss and, to a lesser extent, GIA, both contribute to a sea level fall in Hudson Strait 
(Figures 3d-f), the magnitude of sea level fall is insufficient to compensate for the sea level rise associated 
with Antarctic ice loss. By 2500, the sea level fall across much of the HBC reaches 3 m. Our low-end simula- 
tion suggests a gradual extension of sea level rise into the HBC. By 2500, the extent of sea level rise spreads 
westwards from the Hudson Strait into Hudson Bay, where sea level rise is less than 0.5 m. 
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Figure 4. Total projected sea level change at 2100, 2300, and 2500, relative to 2000. Sum of individual contributions presented in Figure 3. Panels (a-c) 
correspond to the low-end Antarctic ice loss scenario (Golledge et al., 2019). Panels (d-f) correspond to the high-end Antarctic ice loss scenario (Pollard 


et al., 2017). 
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The projection that incorporates the high-end Antarctic scenario is characterized by a transition from sea 
level fall to sea level rise across the HBC over the simulation time. At 2100 under this high-end Antarctic 
scenario (Figure 4d), GIA still dominates in Hudson Bay, James Bay and northern Foxe Basin. In other 
regions of the HBC, the signal of sea level rise from Antarctic ice loss already begins to outpace the signal 
from GIA by 2100. We find that by 2,200 sea level has begun to rise across the entire HBC. Despite the large 
GIA signal in Foxe Basin and James Bay, our projections indicate that these regions would experience up 
to 3 m of sea level rise by 2,200 and that sea level rise in other regions in the HBC reaches up to 4.5 m. In 
2300 (Figure 4e), a maximum sea level rise of 9.5 m is predicted in Ungava Bay in Hudson Strait, with the 
smallest rise of 6.3 m in James Bay. The magnitude of minimum sea level rise under the high-end scenario 
at 2300 is triple the magnitude of maximum sea level fall predicted from the low-end scenario. At 2500 (Fig- 
ure 4f), sea level rise exceeds 14 m in some regions of Hudson Strait. In Foxe Basin, the projected 10 m rise 
by 2500 represents, on average, a 10% increase in water depth, relative to present day. Moreover, the 9 m rise 
predicted for James Bay suggests a 20% increase in water depth relative to present day. As in 2300, we find 
that the magnitude of minimum sea level rise in James Bay and Foxe Basin under the high-end scenario in 
2500 is three times the magnitude of maximum sea level fall calculated from the low-end scenario. 


4.2. Tides 


In Figure 5, we present the M, tidal amplitude changes associated with the total sea level changes shown 
in Figure 4. The M; tidal amplitude responses differ substantially between the low-end (Figures 4a—4c) and 
high-end (Figures 4d and 4e) Antarctic ice loss scenarios. For the low-end scenarios, all three time slices 
show a similar amplitude change pattern with changes increasing in magnitude from 2100 to 2500. Both the 
flood (Figure 5) and no-flood (Figure S8) runs show very similar amplitude change characteristics, therefore 
the description will be based on the simulations permitting flooding. Amplitudes increase in Hudson Strait 
and in the western half of Foxe Basin whereas decreases can be seen along most of the Hudson Bay margins. 


Under the low-end scenario, and in Foxe Basin, where sea level deceases by ~1 m in 2100 and by over 2 min 
2300 and 2500, respectively, occur, the amphidrome in the eastern part of the bay is shifted eastward due to 
increases in the energy fluxes into the eastern part of the bay (see panels a-c in Figures S9 and S10). Small 
decreases in the energy fluxes along the western margin of Hudson Bay can be seen which are associated 
with a small shift of the western amphidrome toward the western margin, thus decreasing amplitudes. The 
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Figure 5. Total projected M, amplitude change at 2100, 2300, and 2500, relative to 2000 for simulations where flooding of low-lying land is permitted with 
increasing sea levels. Panels (a-c) correspond to tidal changes associated with the low-end Antarctic ice loss scenario (Golledge et al., 2019). Panels (d-f) 
correspond to tidal changes associated with the high-end Antarctic ice loss scenario (Pollard et al., 2017). 
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eastern amphidrome shifts north-eastward leading to deceases in amplitudes along the eastern margin. In 
James Bay, the amphidrome in the east moves eastward and becomes degenerate, leading to increases in 
amplitudes in the south-east part of the bay. 


For the simulations adopting the high-end Antarctic ice loss scenario, the changes for the 2100 simulation 
are small and will therefore not be discussed. For the 2300 and 2500 high-end scenarios (as in the low-end 
scenarios), the flood and no-flood responses are very similar. However, in the no-flood simulations the am- 
plitude changes are larger, especially along the eastern margins of Hudson Bay and Foxe Basin (Figure S8). 
In the flood runs, for both 2300 and 2500 (see panels e and f in Figure 5 and S8), large amplitude decreases 
can be seen in Hudson Strait (exceeding 70 cm by 2500), in Ungava Bay and in the southern and northern 
parts of Foxe Basin (exceeding 1 m by 2500 in both regions). In eastern Foxe Basin and along the eastern 
margins of Hudson Bay, amplitudes increase by around 20-30 cm by 2500, whereas along the south-eastern 
margins of Hudson Bay amplitudes decrease by approximately the same amount. 


The amplitude decreases in Hudson Strait for 2300 and 2500 high-end Antarctic ice loss scenarios exceed 
25% and are associated with a strong decrease in energy fluxes through Hudson Strait and into Foxe Basin 
(see panels e and f in Figures S9 and S10). The decrease in energy fluxes together with increases in energy 
losses in the newly flooded areas around the margins of Foxe Basin shift the eastern amphidrome in Foxe 
Basin southward, thus increasing amplitudes in the easternmost part of the basin and decreasing ampli- 
tudes in the south. Similarly, energy fluxes along the western margins of Hudson Bay decrease, shifting 
both amphidromes westward and increasing amplitudes in the east of the bay. In James Bay, the decreased 
energy flux shifts the amphidrome north-westward increasing amplitudes in the south east of the bay and 
leading to decreases in the western half. 


Next, we consider the impact on energy dissipation. For the low-end simulations, integrated dissipation 
over the area of HBC remains similar to present (Table S8). In the no-flood runs, dissipation remains con- 
stant at the present-day value of 232 GW with increases and decreases in Foxe Basin cancelling each other 
out in all three time-slices. In the flood runs, total dissipation increases by 5% in the 2500 time slice. For the 
high-end scenarios, in contrast, the changes are much larger. In the no-flood runs, integrated dissipation 
of the HBC decreases by 23% by 2300 and 40% by 2500. Nearly all of the decreases in dissipation are due 
to reduced energy losses in Foxe Basin, Hudson Strait and Ungava Bay. In the 2500 time-slice, these three 
areas together experience a drop in dissipation of ~80 GW which explains the majority of the dissipation 
decrease. For the flood runs, the dissipation change picture is very similar (see panels e and f in Figures S11 
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and S12), however, the integrated dissipation changes by only 3% and 7% for 2300 and 2500, respectively. In 
these runs, the energy losses in the newly flooded areas nearly balance the dissipation losses away from the 
margins. However, if the newly flooded areas are excluded, the same large dissipation decreases in the three 
basins emerge (see last two columns of Table S8). 


5. Discussion and Conclusions 


We have provided the first regional study of future sea level changes in the Hudson Bay associated with ice 
loss from the Greenland and Antarctic ice sheets and GIA and their impact on ocean tides in the region. We 
showed that both sea level and tides in the HBC are relatively insensitive to Greenland melting but strongly 
sensitive to Antarctic ice sheet evolution, with our simulations indicating a transition from a sea level fall 
to a sea level rise across the HBC that is, dependent on the adopted Antarctic ice loss scenario. In this work, 
we focused on the RCP8.5 emissions pathway, however more moderate warming scenarios such as RCP2.6 
or RCP4.5 could prolong the duration of sea level fall in the HBC. 


Our results indicate that the sea level fall associated with GIA contributes significantly to sea level change 
in the HBC and dominates under low-end Antarctic ice loss scenarios. HBC is dominated by sea level de- 
creases (up to 1.2 and 3 m by 2100 and 2500, respectively). Hudson Strait is the only region that experiences 
sea level rise (up to 0.5 m by 2500). We found that the associated M; tidal amplitude changes are less than 
0.2 m in the HBC by 2500, with tidal amplitude increases in Foxe Basin and western Hudson Strait, and 
tidal amplitude decreases across Hudson Bay, James Bay, and Ungava Bay in Hudson Strait. We have shown 
that small migrations of the amphidromic points are responsible for the observed changes in the M; tides. 


Conversely, GIA does not compensate the projected sea level rise under more rapid Antarctic ice loss in our 
high-end scenario. In the projection adopting the high-end Antarctic ice loss scenario, sea levels begin to 
rise in much of the HBC by 2100, with a widespread transition from sea level fall to sea level rise by 2,200 
and reaching a peak sea level rise of 14.5 m in Ungava Bay in Hudson Strait by 2500. Under the high-end 
scenario, the M; tidal amplitude changes reach over 1 m with the largest decreases occurring in Hudson 
Strait and Foxe Basin and peak increases along the eastern margin of Hudson Bay by 2500. Furthermore, 
as the HBC is presently one of the areas on the planet where the most tidal energy is dissipated (Egbert & 
Ray, 2001), the significant decreases in dissipation in the open ocean areas of Hudson Strait, Foxe Basin and 
Ungava Bay (see panels d-f in Figures S11 and S12) suggested under our high-end scenario simulations 
could affect both local tidal mixing important for biogeochemistry but also the global tidal energy balance. 


Our results suggest that the strong and uniform in sign tidal response of the Hudson Strait, Ungava Bay 
and Foxe Basin system to the sea level variations in both the high- and low-end scenarios could be linked 
to the resonant properties of the basins. A number of previous studies have shown that the HBC straddles 
a number of resonant frequencies close to and either side of the M, frequency (Arbic et al., 2007; Cummins 
et al., 2010; Webb, 2014). Our simulations show that Hudson Strait, Ungava Bay and Foxe Basin change in 
conjunction with each other in response to sea level changes. The moderate decreases in sea level (around 
2 m) seen in the low-end scenario lead to increases in amplitudes, suggesting that the system is moving 
closer to a resonant state. In contrast, the increases in sea level seen in the high-end scenario in the 2300 and 
2500 time slices drive strong amplitude decreases in Hudson Strait, Ungava Bay and Foxe Basin together 
with large decreases in dissipation. In the 2500 high-end simulations (both flood and no-flood), the dissipa- 
tion decreases amount to around 80 GW when newly flooded areas are excluded, nearly halving the energy 
dissipated in the three basins with respect to present day. 


Previous studies investigating resonances have generally considered Hudson Strait and Hudson Bay (to- 
gether with Ungava Bay) as a resonant system (e.g., Cummins et al., 2010). However, our results indicate 
that it is more likely that Foxe Basin and Hudson Strait together with Ungava Bay act as a near-resonant 
system at present; which, perturbed by sea level changes, exhibits large changes in amplitudes and dissipa- 
tion. At present, in Foxe Basin ~110 GW of M, tidal energy are dissipated, whereas in Hudson Bay a mere 
20 GW are lost (Figure S13). Nearly half of the dissipation decreases in the 2500 high-end scenario occur in 
Foxe Basin, and the remainder take place in Hudson Strait and Ungava Bay. This hypothesis is supported by 
the pathway of the energy fluxes into the HBC (Figure S14) and the changes thereof (Figures S9 and S10; 
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see also Figure 1 in Webb et al., 2014) where the main pathway of energy flow into HBC is through Hudson 
Strait and into the eastern part of Foxe Basin with very little energy transfer into Hudson Bay itself. A sim- 
ple back of the envelope calculation of the natural period of Hudson Strait and Foxe Basin (T = 4L / JD J; 
using a distance L of 1,600 km from outside the entrance of Hudson Strait to the south eastern margin of 
Foxe Basin and a mean depth D of 136 m, yields a resonant period T of 48.7 h, which is very close to four- 
times the M, period. A sea level decrease by 2 m as seen in the 2300 and 2500 low-end simulations increases 
the resonant period to 49.0 h, that is, moving it closer to four times the M, period. In contrast, a sea level 
increase of 9 and 13 m as in the 2300 and 2500 high-end simulations results in decreases in the resonant 
period to 47.1 and 46.5 h, respectively, thus moving away from four times the M, period. 


We have carried out additional simulations with the same sea level forcing across HBC but with present-day 
boundary conditions at the open boundaries (not shown) and find very similar results as when we apply 
boundary conditions from the global model forced with sea level changes. A no-flood simulation of the 
high-end 2500 sea level change shows a dissipation decrease of 35% which is very similar to the 40% seen in 
our no-flood runs. We therefore conclude that, to a large extent, the changes seen in Hudson Bay, Ungava 
Bay and Foxe Basin are driven by changes in internal dynamics due the sea level changes rather than by 
external forcing outside our regional model domain. This is in contrast to, for example, the work by Harker 
et al. (2019) carried out for the Australian Shelf. 


Our HBC tidal runs performed at 1/30° degree resolution are more accurate than the global simulations 
shown in Wilmes et al. (2017) with a resolution of 1/8°. They found a global RMS error for M, amplitudes 
of 7.7 cm (global) and 3.8 cm (deep) whereas our new global model setup has reduced errors of 5.6 cm 
(global) and 2.8 cm (deep). However, the responses to sea level changes in their simulations are similar to 
our findings for the HBC simulations. They compared tidal amplitude responses to nonuniform sea level 
changes due to full collapses of the West Antarctic and Greenland ice sheet and a uniform sea level rise of 
12 m. The spatially variable sea level change scenario had a similar magnitude sea level fall over the HBC 
as the 2500 low-end scenario whereas the uniform sea level rise case was most akin to the 2500 high-end 
case shown here. The amplitude responses to the two different scenarios agree in their patterns between 
the two studies which increases our confidence in the results presented in Wilmes et al. (2017), despite the 
different model resolutions. 


Here, we focused on the contributions due to ice mass loss from the polar ice sheets and GIA associated with 
the last deglaciation, however, oceanographic and sea ice changes in the HBC will also impact sea level and 
ocean tides, in particular on timescales shorter than the centennial timescales considered here. Thermal 
expansion and ocean dynamics are expected to contribute on the order of 0.3 m of sea level rise in the HBC 
under RCP8.5 by 2100 (Jevrejeva et al., 2016), similar to the contribution of these effects to globally averaged 
sea level rise of 0.21-0.33 m by 2100 estimated in the IPCC (Church et al., 2013). While regional estimates in 
the HBC are limited on multi-century timescales, the contribution of thermal expansion and ocean dynam- 
ics to globally averaged sea level rise under RCP8.5 is estimated to be between 0.29 and 1.81 m by 2300, and 
0.37-2.77 m by 2500 (Church et al., 2013). We thus expect that the magnitude of oceanographic effects in the 
HBC will be smaller than the contributions from the Antarctic ice sheet and GIA (Figure 3) but could none- 
theless impact the timing of the transition from sea level fall to sea level rise. Incorporating projections of 
variability in ocean dynamics would add to the work presented here, in addition to improving shorter term 
projections of sea level and tide changes in the HBC. Furthermore, we have considered changes in sea level 
and tides in the HBC during the ice-free season. Because Hudson Bay is currently characterized by a long 
sea ice cover season which dampens the magnitude of the tides (Godin, 1986; Kleptsova & Pietrzak, 2018; 
St-Laurent et al., 2008), future work will explore the impact of retreating sea ice on our results. 


Previous simulations of modern tides have had large amplitude errors in the HBC (e.g., Wilmes et al., 2017) 
due to large uncertainty in the bathymetry in the HBC. Because bathymetry in the HBC is poorly con- 
strained, we used tidal simulations and tide gauge data to constrain the bathymetry. We suggested that 
being able to simulate tides more accurately provides an indication of which bathymetry dataset is more 
accurate. We produced a composite bathymetry constrained using TPXO8 solutions and tide gauge obser- 
vations that greatly improves the fit between predicted and observed modern tides, with total, and deep 
ocean root mean square errors around 12 and 3 cm, respectively. Our approach may be used to further 
refine bathymetry models in the HBC and in other regions with sparse data. Our results also highlight the 
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need for high-resolution bathymetry surveys in the HBC, a region that will become increasingly important 
as sea ice cover retreats and access via marine traffic increases in a warming climate. Furthermore, accurate 
bathymetry data has far-reaching implications for residents surrounding the HBC, including, for example, 
future land claim issues related to the formation of land bridges between islands and the mainland (Tsuji 
et al., 2016). 


We emphasize that results from this research have profound impacts on ecosystems and communities sur- 
rounding the HBC. Due to the shallow topography and bathymetry in the HBC, we project that an increase 
in the absolute depth of water of up to 20% could occur by 2500 under our high-end scenario, which would 
lead to changes in the extent of the intertidal zone and shoreline locations, in turn, affecting ecosystems 
including wetlands and marshes. If sea level rise outpaces the rate at which marshes grow, these areas will 
become inundated, reducing the capacity of these ecosystems to attenuate storm surges and flooding events 
(Kirwan & Guntenspergen, 2010). The substantial changes in dissipation and thus tidal mixing have the 
potential to impact water structure and biogeochemistry and thus primary and secondary productivity in 
HBC. Furthermore, the livelihoods, culture, and traditions of indigenous communities surrounding the 
HBC are inextricably linked to the natural environment. Coastal management and policy decisions should 
integrate the physical changes presented here. In doing so, these vulnerable groups will be better equipped 
to deal with future impacts of climate change on coastal regions. 


We conclude that projections of changes in sea level and ocean tides in the HBC should take into account 
not only the sea level changes due to ongoing GIA from the past deglaciation, but also consider a range of 
spatially variable sea level change projections arising from land ice loss, in particular in Antarctica. Uncer- 
tainty in projections of change in coastal areas in the HBC is dominated by uncertainty in projections of the 
evolution of the Antarctic ice sheet. As our understanding of the response of Antarctica to climate change 
improves, so too will our projections of sea level and tide changes in the Arctic. 


Data Availability Statement 


Topography and tidal amplitude data is available through the cited references. Model data is available at 
https://doi.org/10.17605/OSF.IO/UXKFG. 
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